//////////////////////////////////////////////////////////////////////
// MatrixUtil.h
//
// SHEN Fangyang
// me@shenfy.com
//
// Copyright (C) SHEN Fangyang, 2011, All rights reserved.
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
//  Brief Description
//
//	SIMPLE LINEAR ALGEBRA PACKAGE
//
//	Matrix utilities.
//
//////////////////////////////////////////////////////////////////////
#pragma once

#include "SysTypeTraits.h"
#include "boost/numeric/ublas/matrix.hpp"

namespace grflib
{
	namespace LA
	{
		namespace quantify
		{
			namespace ublas = boost::numeric::ublas;

			//quantify matrix to destination range
			template<typename matrix_type1, typename matrix_type2>
			static UINT Quantify(const matrix_type1 &srcMat, 
								 matrix_type2 &dstMat, 
								 typename matrix_type2::value_type threshold,
								 double &scale)
			{
				UINT size1, size2;
				UINT nNonZero = 0;

				//initialize dst matrix
				size1 = srcMat.size1(); size2 = srcMat.size2();
				if (dstMat.size1() != size1 || dstMat.size2() != size2)
					dstMat.resize(size1, size2);

				//find maximum and calculate scale
				double maxval = norm_inf(srcMat);
				scale = w32lib::SysTypeTraits<matrix_type2::value_type>::MaxVal() / maxval;

				//quantify
				for (UINT x = 0; x < size1; x++)
				{
					for (UINT y = 0; y < size2; y++)
					{
						dstMat(x, y) = static_cast<matrix_type2::value_type>((double)(srcMat(x, y)) * scale);
						if (abs(dstMat(x, y)) > threshold)
							nNonZero++;
					}
				}
				return nNonZero;
			}

			template<typename matrix_type>
			static UINT NonZero(const matrix_type &srcMat, typename matrix_type::value_type threshold)
			{
				UINT size1 = srcMat.size1();
				UINT size2 = srcMat.size2();
				UINT nNonZero = 0;

				for (UINT x1 = 0; x1 < size1; x1++)
				{
					for (UINT x2 = 0; x2 < size2; x2++)
					{
						if (abs(srcMat(x1, x2)) > threshold)
							nNonZero++;
					}
				}
				return nNonZero;
			}
		} //namespace quantify

	} //namespace LA
} //namespace grflib